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Abstract 

This paper reviews the theoretical foundation and computational mechanics aspects 
of the recently developed shear-deformation theory, called the Refined Zigzag 
Theory (RZT). The theory is based on a multi-scale formalism in which an 
equivalent single-layer plate theory is refined with a robust set of zigzag local layer 
displacements that are free of the usual deficiencies found in common plate theories 
with zigzag kinematics. In the RZT, first-order shear-deformation plate theory is 
used as the equivalent single-layer plate theory, which represents the overall 
response characteristics. Local piecewise-linear zigzag displacements are used to 
provide corrections to these overall response characteristics that are associated with 
the plate heterogeneity and the relative stiffnesses of the layers. The theory does not 
rely on shear correction factors and is equally accurate for homogeneous, laminated 
composite, and sandwich beams and plates. Regardless of the number of material 
layers, the theory maintains only seven kinematic unknowns that describe the 
membrane, bending, and transverse shear plate-deformation modes. Derived from 
the virtual work principle, RZT is well-suited for developing computationally 
efficient, C°-continuous finite elements; formulations of several RZT-based 
elements are highlighted. The theory and its finite element approximations thus 
provide a unified and reliable computational platform for the analysis and design of 
high-performance load-bearing aerospace structures. 

Keywords: laminated composites, sandwich structures, refined zigzag theory, multi- 
scale, inter-laminar damage, C°-continuous finite elements. 
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1 Introduction 


Lightweight and high-performance characteristics of advanced composite materials 
have spurred a wide range of application of these materials in military and civilian 
aircraft, aerospace vehicles, and naval and civil structures. To realize the full 
potential of composite structures for primary load-bearing components, further 
advances in structural design, analysis methods, and failure prediction and 
progressive damage methodologies are necessary. 

A wide variety of modern civilian and military aircraft use relatively thick 
laminated composite and sandwich laminates for primary load-bearing structures. 
Such structures commonly undergo relatively pronounced design-critical transverse- 
shear deformations and, under certain conditions, thickness-stretch deformations. 
Fail-safe designs require accurate stress-analysis methods, particularly in the regions 
of stress concentration. Moreover, computationally efficient progressive damage 
analysis necessitates accurate modeling of the interlaminar damage modes such as 
delamination. For these reasons, structural theories that account for higher-order 
deformation effects have attracted much attention in recent years. 

Recently, Tessler and co-authors presented an improved structural theory for 
beams and plates, labeled the Refined Zigzag Theory (RZT), that offers substantial 
analytic and computational advantages for the analysis of homogeneous, laminated 
composite, and sandwich laminates [1-5]. This new theory is based on a multi-scale 
approach in which an equivalent single-layer plate theory is used to represent the 
overall {coarse) plate response characteristics and through-the-thickness zigzag 
kinematics are used to model the local (fine ) layer-level behavior. In particular, the 
RZT uses First-order Shear-Deformation Theory (FSDT) as the equivalent single- 
layer theory, and uses sets of piecewise linear continuous functions to model the 
local behavior of the layers comprising a plate. The zigzag kinematic framework 
enables sufficiently accurate and computationally efficient modeling of a wide range 
of homogeneous and heterogeneous laminates without the use of shear correction 
factors. Novel zigzag functions, derived a priori from constitutive relations without 
enforcing debilitating stress-equilibrium constraints, are responsible for overcoming 
several critical shortcomings of the earlier zigzag theories (refer to [1-5] for the 
literature reviews and pertinent discussions). The stress resultants are obtained from 
the equilibrium equations and, as a result, are physically consistent with their 
definitions based on Hooke’s relations. The formulation, which maintains a fixed 
number of kinematic unknowns regardless of the number of material layers, does not 
enforce full continuity of the transverse shear stresses along material-layer 
interfaces, yet is robust. Using the principle of virtual work, equilibrium equations 
and consistent boundary conditions are derived in a variationally consistent manner. 
The variational framework, requiring relatively simple C°-continuous kinematic 
interpolations, provides a convenient means for developing computationally 
efficient and robust finite elements. 

The focus of this paper is to assess the current state of the art in the development 
of RZT and to highlight its recent advances in finite element analysis. To 
accomplish these objectives, this paper first reviews the theoretical foundation of 
RZT for laminated plates and highlights the essential aspects of the piecewise linear 
zigzag functions derived from transverse- shear constitutive relations. Related efforts 
of the zigzag-function enrichment using higher-order polynomials are also 


2 



highlighted. In Section 2, the principle of virtual work is presented, which serves as 
the foundation for developing finite element approximations. Also highlighted in 
this section are (i) a simple and effective method of computing highly accurate 
interlaminar stresses, (ii) the unique modeling of homogeneous plates using the full 
kinematics of RZT, and (iii) a higher-order theory that includes odd and even zigzag 
kinematic terms in the inplane expansions and which also accounts for the thickness- 
stretch deformations. In Section 3, recent efforts to develop RZT-based finite 
element for laminated composite and sandwich beam, plate, and shell structures are 
discussed, focusing on their kinematic approximations, element nodal 
configurations, and modeling capabilities. In Section 4, finite element results for a 
laminated cylindrical shell undergoing large displacements are presented to 
demostrate the latest RZT-based shell modeling capability. Finally, conclusions are 
presented which highlight the salient features of RZT and the finite elements 
developed on its foundation. 

2 Foundation of Refined Zigzag Theory (RZT) 


Consider a multilayer composite plate of uniform thickness 2 h that is composed of 
perfectly bonded orthotropic layers, labelled with superscript (k), where k=l,..., N; 
furthermore, the plate undergoes small- strain deformations under static loading 
while exhibiting negligible inertial effects (refer to Figure 1.) The inplane 
coordinates of the plate are defined by the vector x = (x,,x 2 ) e S m , where S m 
represents the set of points given by the intersection of the plate with the plane z-0 
(the midplane); the symbolism ze[-/i, h] denotes the domain of the through-the- 
thickness coordinate. 



Figure 1. (a) Plate of uniform thickness, 2/z, subjected to transverse pressure 
loading, q , and edge tractions, { T a (a = 1,2) , T_ }; (b) through-thickness notation of 

N -layer laminate. 
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The displacement vector of any material point P(x ] , x 2 , z) in the k th layer is 
defined by the three orthogonal components (u[ k) ,u\ k> ,u_) that are expressed as [4- 

5], 

u\ k} (x, z) = u{x) + z6 x (x) + </\ k} {z.) ¥\(x) 

U { 2 k) (x, z)=v(x) + z0 2 (x) + flf ) (z) ^ 2 (x) 

U z (x,z) = w(x ) ^ 

where u(x) and v(x) are the inplane displacements in the x x and x 2 coordinate 
directions, respectively, and w(x) is the transverse deflection. The symbols <9, (x) 
and 6 2 (x) represent overall bending rotations of a transverse material line element, 
that is initially straight and normal to the midplane, about the positive x 9 and the 
negative x x directions, respectively. These overall rotations represent linearized 
weighted averages of the actual nonlinear through-the-thickness displacement 
variations that are produced when the displacement variations within each layer are 
pronounced. Collectively, these three displacements and two rotations form the basis 
of FSDT and, as such, are viewed herein as the quantities that define the overall, 
first-approximation, coarse displacement characteristics of a plate. 

Refinements to the coarse displacement characteristics, and hence improved 
fidelity, are obtained by introducing additional displacement functions that model 
lamina-scale, or even sub-lamina-scale, responses adequately. In the RZT, the 
refined contributions to the “coarse” inplane displacements are due to the zigzag 

terms (/) {k \z) y/ a {x) {a- 1,2) that appear in Eqs. (1). These refinement functions 
generally produce piecewise-continuous, nonlinear through-the-thickness 
distributions that are defined by the zigzag functions (fi'J ’iz) and their corresponding 
rotation (or amplitude) functions y/ a (x) . The simplest form for <t> f k) (z.) is given by 
a set of piecewise linear, continuous functions with derivatives </> { a k) z that are 

discontinuous at the layer interfaces. Herein, a comma followed by the subscript z 
denotes differentiation with respect to the through-the-thickness coordinate z. This 

particular class of functions is denoted herein by C i0) . 

The localized rotations of the line elements within each layer, which prior to 
deformation are straight and perpendicular to the reference midplane, are obtained 
from Eqs. (1) as 

u al{^z) = d a {x)+/3 (k) y/ a {x) (a =1,2) (2) 

where fi ik) = ^ k \ are constant valued within the k th layer and generally have 
different values for layers with different material properties; thus, they are 
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discontinuous at the layer interfaces. Therefore, the rotations (x, z.) are also 

uniform within the k th layer and generally vary from layer to layer. 

A key step within RZT is to identify its relationship with FSDT. This task is 
accomplished by representing the coarse rotations, 0 a (x ) , as the weighted-average 
rotations, i.e., 

d aW=^\\ u al(*’ z ) dz (« = 1 > 2 ) ( 3 ) 

Substituting Eq. (2) into (3), results in 

\ h _ h Pa ] dz = <f> { a N) (h) - 4 1} i-h) = 0 {a = 1, 2) (4) 

thus guaranteeing equal values of the zigzag functions on the bounding surfaces, i.e., 
(j) f/ N] (h) = ] (—h) (a = 1, 2) . To correlate the displacement of FSDT and RZT at the 

bounding surfaces, z = +h, it becomes immediately apparent that the top and bottom 
values of the zigzag functions must vanish identically, i.e., 

<p { a\h)=$\-h)= 0 (a = 1,2) (5) 

Equations (3) also imply that the following weighted-average transverse shear 
strains are those which correspond to FSDT, i.e., 

/a= ih $-h r ^ dz = w - a + ° a 1 2) (6) 

where henceforth the notation (•) a = d(») / Gx a denotes partial differentiation with 
respect to the midplane coordinate, x a . 

It can now be readily ascertained that within RZT only the coarse kinematic 
variables define the values of the inplane displacements at top and bottom surfaces, 
i.e., 

u\ l) (x,-h) = u(x)-h #j(x), u[ N \x,h ) = u(x) + h <9,(x) 

(x, - h ) = v(x) -h0 2 (x), u[ N> (x, h ) = v(x) + h 0 2 (x) (7) 

With this insight, it becomes clear that the zigzag kinematics can indeed be regarded 
as local ply-level perturbations to the overall, coarse displacement fields. From Eqs. 
(7), the physical interpretations of the coarse variables { u , v , 0 X , 0 2 } in terms of 

the top-surface inplane displacements u at (x) = u ( a N) (x,h) and the bottom-surface 

inplane displacements -h) also become apparent, i.e., 
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( 8 ) 


w(x) = - [ lly (x) + u lb (x)~| , e x (x) = - 1 - 1 " Uy (x) - u ]b (x) 

2 L J 2 li L J 

v(x) = -\u 2t (x) + u 2b (x)], 0 2 (x) = — \u 2t {x)-u 2b {x)\ (9) 

2 L J 2/7 L J 

These same relations hold for these kinematic variables used in FSDT. 

For certain classes of problems, modelling advantages can be exploited by using 
a different set of primary unknown functions instead of the seven appearing in Eqs. 
(1). For example, in finite element analyses that use the sub-laminate concept [6-9], 
the inplane displacements of the top and bottom bounding surfaces of a sub-laminate 
are particularly useful as primary unknowns. Thus, by treating a laminate as a 
collection of contiguous sub-laminates, the RZT equations presented herein can be 
modified as follows. From Eqs. (8) and (9), it is seen that an alternate description of 
the RZT displacement fields is obtained by expressing the four coarse variables { u , 
v , 0 X ,6 2 } in terms of the four inplane displacements of the top and bottom surfaces 

{u at (x), u ab (x) } (a = 1,2), and then substituting the result into Eqs. (1). In this 
way, RZT is redefined in terms of the alternative set of seven displacement 
variables: { u at (x ), u ab {x ) , y/ a {x), w(x ) } (a = 1,2) which can be applied on a 

laminate and a sub-laminate level. Because this transformation of the primary 
unknowns is well defined, this version of RZT is expected to have the same analytic 
accuracy as the original RZT. 

Using the linear strain-displacement relations of elasticity theory and the 
displacement assumptions given by Eqs. (1), the RZT strains are given as 

£ u = u ,\ + zO n +(f{ k) Y u 
s 22 = v 2 + z0 22 + (jh Yl2 

Y\2 = U 2 + V| + Z.{0\2 + ^2,]) + ^ ,< V],2 + ^ V 2 J (10) 

Yal =Ya+Pa ] Ya (« =1 > 2 ) 

where it is noted that the inplane strains, { e[ k x ,s 22 , /l 2 }, are piecewise linear 

through the laminate thickness. The transverse shear strains, y ( a k '} , are piecewise 

constant , i.e., they are constant within each material layer but discontinuous along 
the layer interfaces, in contrast to those of FSDT which are constant across the total 
laminate thickness. 

Each layer of the plates considered herein is presumed to be specially orthotropic 
with respect to a set of principal material coordinate directions. In addition, each set 
of layer principal material axes have two axes that reside in the plane of the layer 
and a third that is perpendicular to that plane. This third axis is parallel to the z-axis 
for the plate coordinates defined herein. The two principal axes that reside in the 
plane of the layer are generally noncoincident with the inplane coordinate axes of 
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the plate. In terms of the plate coordinate system, the generalized Hooke's relations 
for the k th layer are given herein by 


^11 

(■ k ) 

1 

'sD 

(N 

1 

(k) 


Cr 22 

> = 

^22 ^26 

< 

8 22 

r 12 


sym. C 66 


712, 


(k) 


or o (t) =C (t) E (t) 


for the inplane stresses and by 


i r 


'2 z 


n z. 


Q22 

sym. 


Qn 

Qn 



or 


T (k) = q(A-) (k) 


( 11 . 1 ) 


(11.2) 


for the transverse shear stresses. The symbols C\ k) and Q - k ' denote the transformed 

elastic-stiffness coefficients referred to the plate coordinate system. In RZT, these 
coefficients are also based on the presumption that the transverse-normal stresses are 
negligible. 

The next step in formulating RZT is to provide a mathematical description of the 
zigzag functions. Thus, for the k lh material layer, which is located in the range 
[z (k 1 ,,%, ]- the zigzag functions are given as (see Figure 1(b) depicting the 
lamination notation), 

+ + <a = 1.2) (12.1) 


where 


£■“’=[(*- z a - n )/lt lk) - 1] s [-1, 1] (k = 1,..., N) 


( 12 . 2 ) 


with the first layer, k - 1 , beginning at z {(}) = -h , the last (V th layer, k = N , ending 
at z. (N} = h , and the k th layer ending at z (k) = z (k _ u + 1h (k) , where 2 h (k) denotes the 
k th layer thickness. Evaluating Eqs. (12) at the bottom (C ik> = -1) and top (C ik) = 1) 
surfaces of the k* layer, gives rise to the definitions of the interface displacements 

=4 t) (4' m =-l). 4«)=4‘ ) ((T , ‘’=l) (k = l,...,N) (13) 

whereas according to Eqs. (5) the zigzag functions vanish at the bottom 
C k = 1, £ (]) = -1) and top ( k = N, £ (N) = 1) plate surfaces, i.e., (j) aiQ) = (p a(N) = 0 . 
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(b) Zigzag functions (j. \ (k \z ) and (/h'iz) . 


Figure 2. Notation for a three-layer laminate and </> ( a k) (z) zigzag functions defined 
in terms of interface displacements, (j) a{k) . 

The remaining interface displacements can be obtained from the simple expression 

< i4 > 

i=l 


In [5], two alternative expressions have been derived for J3^ ( k =1,..., TV — 1 ); they 
are given as 


(a) Pa = 


1 


y* k) 

h hQ ( , k) 


-1, (b) A 


(k) = . 

a 


Ak) 




-1 (a =1,2) (15) 


k= 1 
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where and S ( a k J ( a = 1, 2) denote, respectively, the diagonal transverse-shear 

stiffness and compliance coefficients of the k th layer, respectively. Note that these 
expressions are only slightly different from each other for the general lay-up case, in 
which the transverse-shear constitutive matrix has a non-zero off-diagonal term ( 

Q [ 2 ^ 0). The two definitions, however, are identical for the special, decoupled case 

for which Q { 12 } = 0. From the predictive perspective, either form of (3 (k) yields 
excellent results, with the second form exhibiting slightly improved results, 
particularly for sandwich plates. 

Further accuracy enhancements within the present RZT methodology can be 
achieved by using higher-degree polynomial expansions for the zigzag functions 

<j > (k) , that are piecewise continuous nonlinear functions, while keeping the basic 
kinematic assumptions, Eqs. (1), unchanged (e.g., refer to [10-13]). 

2.1 Variational framework 


The principle of virtual work, for the case of negligible body forces and zero shear 
tractions on the top and bottom bounding plate surfaces, may be written as 

h 


IK a (k)T S£ (k) +T k)T Sj (k> 'jdz. clS - J qdwdS 

S m -h S m 


n 

■JJ[ 


+ T 2 Su [ 2 > +T z Su z 


,(*) 


ds dz= 0 


(16) 


C„~h 


where 5 is the variational operator; C 1 is the applied transverse pressure attributed to 
the middle reference surface, S m ; [ T a (a = 1,2) , f_ | are the inplane and transverse 
shear tractions that are prescribed along the cylindrical edge surface of the laminate, 
S a =C a xs . The bold-face vectors appearing in equation (16) are defined by Eqs. 
( 11 ). 

Integrating Eq. (16) across the laminate thickness, while accounting for the 
relationships in Eqs. (1) and (10), yields the corresponding two-dimensional 
statement of the principle of virtual work 


[ [N’ / „d'e m + M T b Se b + Q'd'e v - qSw] dS - [ F c)u ds = 0 

J S m ‘ ‘ J 


(17.1) 


where the membrane stress resultants and conjugate strain measures are defined as 
K={N 1 ,N 2 ,N l 2 } = ljatf,<T£\T™}dz, *1 ={u A ,v 2 ,u 2 +v x } (17.2) 
Likewise, the bending stress resultants and conjugate strain measures are defined as 


9 



(17.3) 


Mfo = |Mj, Mf , 


m 2 ,m!,m 12 ,m ? 2 , 



= J_* ^l (A) ^ifi o-22 } > zr[ 2 \(t)[ k) T$\ 4 k) h [i}dz 

e b = {^1,1’ V0,1’ ^2,2’ ^2,2’ ^1,2 + ^2,1’ ^1,2’ V^2,l } 


The transverse shear stress resultants and conjugate strain measures are defined as 


Ql - {g 2 , Qt, Go Qt } = /g fc) rg\ A (A) *£>}& 

e I ={ M; > 2 + ^2^2> W ,\ +0 v¥x} 


(17.4) 


The force and moment resultants due to the prescribed tractions and their conjugate 
displacements are defined as 


F r - { N ln , N 2n , Q zn ,M h „ M ln , Mi , M*, } 

= \ h _ h {T V T 2 ,T Z , zT v zf 2 , (tl k %,(f^f 2 }dz (17.5) 

u T = {u,v,w,0 l ,0 2 ,i// l ,i// 2 } 

All elements of Eqs. (17) with the superscript </> are associated with the zigzag 
functions. 

The stress resultants are readily obtained in terms of the two-dimensional plate 
strain measures by integrating the expressions in Eqs. (17) through the laminate 
thickness, resulting in the following constitutive relations of RZT: 


n; 


A B O' 

e m 

M. 

.= 

B 1 D 0 

‘ e ;. 

Q SJ 


0 0 G 

, e s. 


(18) 


where A = [A (/ ] 3x3 denotes the membrane stiffness matrix, B = [ B- ] 3x7 is the 
membrane-to-bending coupling stiffness matrix, D = [D (/ ] 7x7 is the bending stiffness 
matrix, and G = [G ;/ ] 4x4 is the transverse shear stiffness matrix. The explicit forms 
of these stiffness coefficients are found in [4,5]. Note that for a general composite 
lay-up, the membrane-to-bending coupling stiffnesses are nonzero; that is, B jj ^ 0 . 
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2.1.1 Special case 


When the zigzag functions in the x, and x 2 directions are identical, 

(f>[ k \z) = <f> 2 ) (z) = (/) {k \z) (k N) . As a result, the two twisting moments 

associated with the zigzag kinematics are given by the same expression (refer to Eqs. 
(17.2) and [14]) as follows 



This simplification gives rise to the following reduced form of Eqs. (17.2) 

M T h ={M l ,M 2 ,M l2 ,Mf ,M%,Mf 2 } 



— {^ 1 , 1 5 ^ 2,2 ’^ 1,2 + ^ 2,1 > 1 ^ 1 , 1 > 1 ^ 2,2 > 1 ^ 1,2 + ^ 2 . l ) 

and the resulting simplifications in Eqs. (18). These constitutive relations provide a 
clear physical interpretation of the new quantities associated with the zigzag 

kinematics. Specifically, {Mf,M%,Mf 2 } represent the bending and twisting 
moments due to the zigzag related cross-sectional distortions of the layers, whereas 
{^ 11 ,^ 22 ^ 12 +^ 21 } are their conjugate curvatures due to the normal zigzag 

rotations, i//, and t// 2 . 

2.2 Interlaminar stresses 

An important attribute of any plate theory that is formulated to be robust at the local 
level is the ability to predict adequately interlaminar stresses that can be used to 
assess laminate failures. To provide adequate estimates of interlaminar stresses, it is 
often customary to integrate the three-dimensional equilibrium equations of 
elasticity theory. In this approach, the interlaminar transverse shear stresses are 
obtained by integrating the jq and x 2 derivatives of the inplane stresses 




Z ( (7 a '> 4. T ^ t 
V°in + M2, 2' 



+ T (k) 
22,2 + M2 


1 ) 


dz 


( 21 . 1 ) 


whereas the interlaminar transverse normal stresses are obtained by integrating the 
derivatives of the transverse shear stresses as 


z 



-h 


+ T 


(*) 

2z,2 


)dz 


(21.2) 
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Although this post-processing procedure has been widely advocated for use with 
FSDT and many higher-order theories, the adequacy of this scheme can be seen to 
be directly linked to (a) the accuracy of the inplane stresses obtained from 
constitutive relations of the underlying theory, and (b) the accuracy of the 
derivatives of these stresses (or second derivatives of the kinematic variables, refer 
to Eq. (10)), with this issue particularly relevant when finite element analyses are 
performed. 

Regarding issue (a), it has been observed that FSDT and the vast majority of 
available structural theories often underestimate the inplane stresses; this is 
especially the case when highly heterogeneous and relatively thick laminates are 
analyzed, including the sandwich construction. Relative to issue (b), reasonably 
good improvements in estimating derivatives of the inplane stresses can be achieved 
by using smoothing techniques, e.g., [15]. 

Application of Eqs. (21), using RZT inplane stresses, has been shown to be 
highly effective, yielding accurate interlaminar stresses that are comparable to those 
predicted by three-dimensional elasticity. This robustness exists because the 
piecewise linear inplane stresses of RZT are consistently accurate, even for highly 
heterogeneous and relatively thick composite and sandwich laminates [1-5]. 

2.3 RZT modelling of homogeneous plates 

The key property of a zigzag function is that it vanishes identically when the 
material has homogenous transverse shear properties across the total thickness. 
Consequently, when the zigzag terms vanish identically, the kinematic assumptions 
given by Eqs. (1) revert to those of FSDT, which require transverse shear correction 
factors. Tessler et al. [5] have found a simple and effective way to use the full 
power of RZT’s kinematic field to model homogenous plate (and beam) problems, 
without the need for transverse shear correction factors. In this approach, a 
homogeneous plate is modelled as a multilayer heterogeneous laminate, in which the 

transverse shear moduli, G {k 3 (a = 1, 2) , vary only slightly from layer to layer, i.e., 

Ga3=G a3 (\ + £ (k> ) (« =1,2; k =l,...,N) (22.1) 


where (a = 1, 2) denote the constant values of the shear moduli, and where 

6-1 1 = + Z (k) + Z (k-l) Z {k)] ^ 1 (s^l) (22.2) 

is a layerwise, dimensionless coefficient that is a function of the position of the k th 
layer within the laminate thickness. Thus, since s (k) <§C 1 (e.g., by setting s = 10 4 ), 
the values of G [k 3 are only slightly perturbed compared to the constant shear- 
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modulus values C a3 of the corresponding homogeneous plate. Thus, in this 
approach, the homogeneous plate is replaced with a corresponding plate with an 
infinitesimal degree of heterogeneity, which for all practical purposes is 
homogeneous. This approach is effective even when the deviation from the constant 
value of the shear modulus of the homogeneous plate is less than 1/100 of a percent. 
The results in [5] have shown that as the number of material layers increases within 
this infinitesimally heterogeneous laminate, the solution approaches that of a 
homogeneous plate, exhibiting the well-known parabolic distribution of the 
transverse shear stresses across the plate thickness, as well as the cubic distribution 
of the inplane stresses, with the latter being particularly evident in relatively thick 
plates. Remarkably, the solution itself finds the correct transverse shear stress 
distributions, albeit in a piecewise constant manner, thus validating the lack of 
necessity for the transverse shear correction factors that are commonly used within 
FSDT and even within some theories of higher order. 

2.4 Higher-order RZT including transverse normal deformations 

The RZT presented herein incorporates FSDT as the underlying baseline plate 
theory that models the overall, coarse plate behavior. Thus, it follows that higher 
order RZT theories can be formulated by using other underlying baseline plate 
theories with the zigzag kinematics used herein or with zigzag kinematics of higher 
fidelity. For example, Barut et al. [14] combined an underlying { l,2}-order baseline 
theory, that uses an overall parabolic transverse displacement assumption, with 
piecewise quadratic RZT-based zigzag inplane displacements to obtain a higher- 
order refined zigzag plate with eleven kinematic variables; four more than the RZT 
presented herein. The displacement assumptions for this theory, labeled as RZT 12 2) , 
are given as 


u[ k ) (x, z) = u(x) + Z0 X (x) + ’ ( z ) W\ 1 ( x ) + (z) Vn^) 

h 

U { 2 k) (x,z) = v(x) + Z 0 2 (x) + </>¥ ) (z)l// 2 l (x) + j</>2 ) (z)l// 22 (x ) 

h 


(*)/ 




(23) 


u (x,z) = w(x) + - Wj(x) + 
h 


f z 2 O 


w 2 (x) 


where (f) ik) are the piecewise linear zigzag functions within the k th layer, adopted 
from the original RZT formulation presented herein; y/ aX (x) and y/ a2 (x) (a = 1,2) 
are the four zigzag amplitudes that permit both non-symmetric, y/ al (x), and 
symmetric, y/ a2 (x), inplane zigzag deformation modes. In addition, an average 
transverse normal stress, cr_, , is independently assumed as a cubic function through 
the laminate thickness, in the form proposed by Tessler [16]; that is 
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(24.1) 


o\ 7 (x, z) = CT-o (x) + cr 7l (x) 


f 3 A 

z z 


'zl' 


h 3 h 3 


This representation of the normal through-the-thickness stress yields a parabolic 
expression for its derivative cr__ ,(x, z) , given by 


cr 77 _ 7 (x,z) 


_ q~ z i( x ) 
h 


V /?2 y 


C7__(x,Z = ±/l) = 0 


(24.2) 

(24.3) 


Equations (24.3) constitute the exact conditions on cr__ in the sense of three- 
dimensional equilibrium equations of elasticity theory for the case of zero-valued 
transverse-shear tractions on the bounding surfaces. The cr_ 0 (x) and cr_, (x) terms in 

Eq. (24) are expressed in terms of kinematic variables of the theory by way of the 
least- square transverse strain compatibility relations originally proposed in [16]. 

Thus, in this new theory, both transverse shear and transverse normal (thickness- 
stretch) deformations are included. The present theory models accurately the in- 
plane and transverse stress components through the thickness. The addition of the 
even (symmetric) inplane zigzag modes, which are not included in the seven- 
variable RZT, contributes to some improvements of the inplane displacement, strain 
and stress response, with further improvements in the interlaminar transverse shear 
stresses. Furthermore, RZT 12 2 ’ appears to be an excellent candidate for developing 
efficient and accurate C°-continuous finite elements. 

3 Finite element approximations 


The variational statement given by Eq. (17.1) can be integrated by parts to produce a 
consistent set of equilibrium equations and boundary conditions, resulting in seven 
partial differential equilibrium equations in terms of seven kinematic variables [4-5]. 
The equilibrium equations can be solved exactly or approximately depending on the 
complexity of the material lay-up, boundary conditions, and loading. Refer to [1-5] 
for details of analytic and approximate solutions of simply supported and 
cantilevered beams and plates made of laminated composite and sandwich 
construction. 

Alternatively, to enable large-scale analysis of complex aerospace structures, the 
variational statement, Eq. (17.1), can be discretized using finite elements. As in 
FSDT, the strain measures within RZT are given as first-order partial derivatives. 
The direct implication is that computationally efficient C°-continuous beam, plate, 
and shell finite elements can be developed, and that the legacy finite element 
technology that has been developed for FSDT can also be used with RZT. 
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3.1 RZT beam elements 


Recently, Gherlone et al. [17,18] and Onate et al. [19] derived planar beam finite 
elements based on the RZT beam formulation in Tessler et al. [1]. When the beam 
deformations are restricted only to the x l - z plane, the second equation in Eqs. (1) is 

identically zero, and the theory reduces to the two displacement components, n\ k) 
and u z . Consequently, there remain four independent kinematic variables, 
u(x l ), wiv, ), 6 X (x t ), and i// l (v, ) , that describe the membrane, bending, transverse 
shear, and zigzag deformations. 

Gherlone et al. [17,18] derived several low-order RZT beam finite elements with 
the aim of achieving the best compromise between accuracy and computational 
efficiency. The four kinematic variables use the anisoparametric (aka 
interdependent) interpolations, where the polynomial degree of w(v, ) is one order 

higher than those approximating the u{ jq), <9, t.q) , and y/\(x { ) variables. Such 
interpolation strategy enables free of shear locking element behaviour for slender 
beams. With an initial assumption of a parabolic distribution for w(x ] ) , the authors 
explored several shear-related constraint conditions that gave rise to two-node 
elements and a coupled form for the w{x x ) interpolation. The constraint condition 
requiring a constant variation of the transverse shear force gave rise to a remarkably 
accurate two-node beam element. For further details of this formulation and for 
elastostatic solutions of simply supported and cantilevered beams of various 
slenderness ratios and lamination properties, the reader is referred to [18]. 

Onate et al. [19] derived a two-node, RZT beam finite element using linear, 
isoparametric interpolations for the four kinematic variables, u{x x ) , vv( .v, ) , 0 X {x x ) , 
and y/ x (v, ) . To achieve an element capable of modelling slender beams without 

incurring the shear locking effect, the authors used reduced integration (a one point 
Gaussian quadrature rule) of the transverse-shear strain energy of the beam element. 
One of the most interesting numerical solutions presented was for a laminated beam 
with an imbedded delamination that was modeled as a compliant layer. It was shown 
that the RZT beam element is computationally efficient and is capable of high- 
fidelity modelling of imbedded delaminations. 

3.2 RZT plate and shell elements 

Versino et al. [20,21] developed six- and three-node triangular RZT -based plate 
finite elements. Adopting linear shape functions for the in -plane displacements, 
bending rotations, and zigzag amplitudes, and a quadratic shape function for the 
transverse deflection, (i.e., using the Tessler-Hughes anisoparametric interpolation 
strategy [22]), the element interpolations in terms of the linear area-parametric 
coordinates L ( have the form 
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3 3 

u(x l , x 2 ) = Y j U jLj , v(xj ,x 2 ) = X VjLj , w(x l , x 2 ) 

i= 1 i=l 

3 3 

6>j C*j, x 2 ) = ^ °\ i L i - 6' 2 (x,,x 2 ) = ^ d 2l L t 

i= 1 r'=l 

3 3 

Ul <*1 > *2 ) = Z A' ’ ^2 (*, > *2 ) = Z ^2 i L i 

i= 1 i=l 

where i e {1,2,3} is an index ranging over the three comer nodes; 
k e {l,m l2 ,2,m 23 ,3,^3i } ranges over the corner and mid-edge nodes; and P k defines 
the standard set of quadratic shape functions (see the nodal configurations in Figure 
3.) The element topology includes six nodes; however, the mid-edge nodes have 
only degrees-of-freedom (dof) associated with the transverse deflection. Using these 
interpolations, the six-node element, called unconstrained, is readily derived by 
introducing Eqs. (25) in the variational principle, Eq. (17), while carrying out the 
necessary matrix and variational operations to obtain the element stiffness equations. 

In addition, the authors derived a three-node, constrained anisoparametric 
element whose interpolation functions for the deflection variable in Eqs. (25) are 
constrained by suitable edge constraints, analogous to those which have been 
explored by Tessler and Hughes [22] for plates modeled with FSDT, and by 
Gherlone et al. [17,18] for beams modeled with RZT. The explicit edge-constraint 
procedure enforces the mid-edge w dof to be dependent on the comer-node dof of 
the rotation variables of the corresponding edges. By imposing one-dimensional 
constraint relations, requiring the two shear strain measures to be constant along the 
element edges, there results a coupled-form, parabolic deflection that is compatible 
along the element edges given by [21] 

3 3 

w(x { ,U) = X L < h ', + X[( 0 u + Wu ) L U+{ °2< + CWli ) L 2i ] 

i= l ;=i 

with 

A; = ^ (bk Lj ~bjL k ), L 2i = UjL k - a k L- ) 

a i= X lk- X lj > b i= X 2j~ X 2k 

where the subscripts are given by the cyclic permutation of i e [1,2,3] , j e [2,3,1] , 
and k e [3, 1,2} ; and where c is either 0 or -1, depending on the constraint strategy 
used in the element formulation. 


( 26 . 1 ) 

( 26 . 2 ) 


=1 


w, 


k P k 


k =i 


( 25 ) 
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Figure 3. Six-node (unconstrained) and three-node (constrained) RZT -based 
anisoparametric plate elements. 

An extensive numerical study was carried out in [21] on various laminated 
composite and sandwich plates undergoing elasto-static deformations. The results of 
this study demonstrate that the RZT-based anisoparametric elements possess 
superior accuracy and excellent convergence characteristics over a wide range of the 
span-to-thickness ratio. In addition, it was shown that the three-node elements 
provide the best compromise between computational efficiency and accuracy. 
Furthermore, the RZT models provide superior through-the-thickness predictions of 
the inplane and interlaminar stresses over comparable FSDT models, especially for 
relatively thick and highly heterogeneous laminated plate and sandwich plate 
constructions. 

The latest finite element implementation of RZT, by Versino [23], is focused on 
developing robust anisoparametric shell finite elements that include the drilling dof 
(i.e., the dof associated with a rotation about the normal to the element’s planar 
surface) and geometrically nonlinear deformations. In the next section, the 
application of RZT to an elastic cylindrical shell undergoing large displacements 
under quasi-static loading is demonstrated. 

4 Nonlinear deformations of a laminated cylindrical shell 

In this demonstration problem, a cylindrical shell is subjected to a quasi-static 
concentrated force, F, applied at point C (Figure 4). The global shell dimensions are 
L=254 mm, r=2540 mm, and 9=0. 1 rad. The shell wall is a three-layer laminate with 
a total thickness 2h=24 mm and with each layer having the same thickness. The 
prescribed boundary conditions are such that the straight edges are fully clamped 
and the curved edges are free. The mechanical material properties of isotropic 
material layers A and B are summarized in Table 1, and the laminate stacking 
sequence is given by [A/B/A]. Note that the elastic modulus of the middle layer is 
two orders of magnitude less than those of the top and bottom layers. This aspect, 
typical of sandwich construction, presents a significant challenge for any structural 
theory, since the through-the-thickness distributions of the inplane displacements 
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exhibit interface slope discontinuities that are not accounted for in most structural 
theories. 

To establish a viable reference solution for this problem, a high-fidelity, three- 
dimensional geometrically nonlinear FEM analysis was performed using the 
ABAQUS commercial code [24]. Depicted in Figure 4(a) is an FEM model which is 
arrived at from a convergence study. The model is based on C3D20 brick elements, 
discretizing a symmetric quadrant of the shell, using a 64x64x6 mesh (64 elements 
along each edge and 2 elements across the thickness of each layer). In addition, two 
shell-based geometrically nonlinear solutions were obtained: (a) an ABAQUS 
solution using the S3 (three-node, FSDT) shell elements, and (b) an RZT-based 
solution using RZT3 - a three-node anisoparametric element with drilling dof. Both 
shell models are based on the 32x32 mesh subdivisions that span the shell’s 
symmetric quadrant, as shown in Figure 4(b). 

Comparisons of the three nonlinear FEM solutions are shown in terms the load- 
deflection curve at point C in Figure 5. It is evident from these results that the three- 
dimensional FEM solution and the corresponding RZT3 solution are in very close 
agreement over the entire range of the applied loading. In contrast, the S3-model 
predictions are considerable less accurate, particularly over the range of larger 
displacements. 

Table 1. Mechanical material properties for isotropic layers. 


Material 

Young’s 
modulus, 
E (MPa) 

Poisson 

ratio, 

V 

A 

31 x 10 2 

0.3 

B 

31 



(a) Solid model (b) Shell model 

Figure 4. Finite element discretization for cylindrical shell. 
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Figure 5. Applied force vs. deflection at point C. 


5 Conclusions 

In this paper the latest advances in the development of a shear-deformable, Refined 
Zigzag Theory (RZT) for homogeneous, laminated composite, and sandwich beams 
and plates have been reviewed and the salient features highlighted. This theory 
maintains four kinematic unknowns for beams and seven for plates, regardless of the 
number of layers, material composition, or layer stacking sequence. In addition, it 
has been shown that the theory can be expressed in terms of displacement quantities 
that are amenable to sub-laminate modelling techniques. A higher-order plate theory 
that incorporates the RZT presented herein has also been described, which includes 
the thickness-stretch deformations in addition to the transverse shear deformations, 
and has eleven kinematic variables. This higher-order theory illustrates how the 
basic RZT presented herein can be adapted to obtain results with a very high degree 
of fidelity. 

Analytic and computational aspects of the RZT and several new finite elements 
for beams, plates, and shells have also been highlighted. The unique characteristics 
of RZT and its finite elements is that they permit efficient modelling of a wide range 
of problems including homogenous beams and plates, composite laminates, 
sandwich construction, and compliant-layer delamination modelling. Such solutions 
do not require shear correction factors or any increase in model complexity or 
computational effort. In addition, RZT -based finite elements are amenable to much 
of the legacy finite element technology, especially that for FSDT-based elements. 
Because of these attributes, RZT has considerable promise for becoming the theory 
of choice for many practical applications including the computationally challenging 
problems of progressive damage modelling in composite structures. 
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